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Abstract 

On the basis of quantum Monte Carlo (QMC) simulations we study the formation 
of Mott domains in the one-dimensional Hubbard model with an additional con- 
fining potential. We find evidences of quantum critical behavior at the boundaries 
of the Mott-insulating regions. A local compressibility defined to characterize the 
local phases exhibits a non-trivial critical exponent on entering the Mott-insulating 
domains. Both the local compressibility and the variance of the local density show 
universality with respect to the confining potential. We also study the momentum 
distribution function of the trapped system, and determine its phase diagram. 

Key words: Optical lattices, metal-Mott-insulator transition, Degenerated Fermi 
gases, Strongly correlated systems 
PACS: 03.75.Ss, 05.30.Fk, 71.30.+h 



1 Introduction 

The study of ultracold quantum gases has become in the last decade a field of 
intense experimental and theoretical research [1,2,3,4]. Bose-Einstein conden- 
sation (BEC) was the main motivation starting the intensive analysis of such 
systems. It was observed for the first time in a series of experiments on dilute 
vapors of alkali atoms cooled down to extremely low temperatures (fractions 
of microkelvins) [5,6,7]. 

Recently, new features that allow to go beyond the weakly interacting regime 
and access strongly correlated limits have been added to the experiments. 
Particularly important, due to its relevance for condensed matter physics, 
has been the introduction of optical lattices. They allowed the study of the 
superfluid-Mott-insulator phase transition in three-dimensional [8], and one- 
dimensional (ID) [9] systems. The presence of the optical lattice and the fact 
that the particles interact only via contact interaction, lead in a natural way 
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to the Hubbard model as a paradigm for these systems. A theoretical work 
proposing such experiments [10], and quantum Monte Carlo (QMC) simula- 
tions [11,12,13] have examined these systems in detail. It has been found that 
incompressible Mott insulating phases appear for wide ranges of fillings, and 
always coexist with compressible phases, so that local order parameters [11,13] 
have to be defined to characterize the system. 

In the fermionic case, the metal-Mott-insulator transition (MMIT) [14,15] 
has not been observed yet. However, recent experiments succeeded in loading 
single species ultracold fermionic gases on an optical lattice [16,17], allowing 
the realization of an ideal Fermi gas on a lattice. Progress in this field, loading 
more than one component fermions and reducing the occupation per lattice 
site, could lead to the realization of the metal-Mott-insulator transition on 
optical lattices. 

Motivated by this expectation we studied, using QMC simulations, the ground 
state of the one dimensional (ID) fermionic Hubbard model with a harmonic 
trap [14,15]. In the present work we review those references adding some new 
results and improving others. As in the bosonic case [11,12,13], we find that 
Mott domains appear over a continuous range of fillings and always coexist 
with compressible phases. Therefore, we define a local order parameter, that 
we call local compressibility, in order to characterize the local phases present 
in the system. By means of this local compressibility, we analyze in detail the 
interface between the metallic and insulating region finding that critical be- 
havior sets in, revealing a new critical exponent. Furthermore, the behavior of 
the local compressibility and the variance of the density are found to be uni- 
versal in this case independently of the confining potential and the strength 
of the interaction. Hence, universality appears as usual in critical phenom- 
ena. The momentum distribution function, a quantity usually accessed in the 
experiments, is studied in detail. The results obtained show that due to the in- 
homogeneous character of these systems, this quantity does not exhibit a clear 
signature of the formation of Mott domains. Finally, we obtain the phase dia- 
gram for fermions confined in ID harmonic traps. It allows to compare systems 
with different fillings and curvatures of the confining potential, so that it could 
help to understand future experimental results. 

The exposition is organized as follows. In the next section, we give a short 
introduction to the Hubbard model on periodic systems. In Sec. 3, we study 
local quantities within the Hubbard model with an additional harmonic trap. 
In particular, we define a local compressibility that acts as a local order param- 
eter to characterize the Mott-insulating phases. In Sec. 4, we study in detail 
the region where the transition between the metallic and the Mott-insulating 
phases occurs. The momentum distribution function is analyzed in Sec. 5. In 
Sec. 6, we discuss the phase diagram for these systems. Finally, the conclusions 
are given in Sec. 7. 
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2 The Hubbard model 



The Hubbard model in a periodic lattice [18,19] has been intensively studied 
in condensed matter physics as a prototype for the theoretical understand- 
ing of the MMIT. This model considers electrons in a single band, and its 
Hamiltonian can be written as 



H 



C icr C i+ltj 



H.C.J + Uj2 n it n il 



where c- CT , c la 



,t 



are creation and annihilation operators, respectively, for elec- 
trons with spin a at site i, and n icr = c\ a c ia is the particle number operator. 
Electrons are considered to hop only between nearest neighbors with a hopping 
amplitude t, and the on-site interaction parameter is denoted by U (U > 0). 
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Fig. 1. General schematic phase diagram for the metal-Mott-insulator transition 
obtained from the Hubbard model. The shaded area represents a metallic region 
that is under the strong influence of the MMIT. The routes for the MMIT (see 
text) has been signaled with arrows. 

It is well known that the Hubbard model [Eq. (1)] displays MMIT phase 
transitions at half filling and at finite values of the on-site repulsive interaction 
U [20]. (In the case of perfect nesting 1 the transition occurs at U = 0.) 
The schematic phase diagram is shown in Fig. 1, in the plane U/t vs n. The 
two routes for the MMIT are the filling- controlled MMIT (FC-MMIT), and 
the bandwidth-controlled MMIT (BC-MMIT) [20]. Metallic regions very close 
to the Mott-insulating phase [shaded area in Fig. 1] are under the strong 



1 Nesting refers to the existence of parallel sections in the Fermi surface. This means 
that excitations of vanishing energy are possible at finite momentum. Perfect nesting 
appears for periodic systems in ID at any filling since in ID the Fermi surface 
consists only of two points. In addition, it appears for hypercubic lattices at half 
filling since in those cases the Fermi surface is a hypercube. 
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influence of the MMIT, and in general exhibit anomalous features like mass 
enhancement and carrier number reduction [20]. In ID systems, the metallic 
phase (Luttinger liquid) is characterized by gap less excitations and by a 2fc f 
(where kp is the Fermi momentum) modulation of the spin-spin correlations. 
On the other side, in the Mott-insulating phase (n = 1) the system exhibits a 
charge gap, and quasi-long range antiferromagnetic correlations. At n = 0, 2 
the system is in a trivial band insulating state. 



3 The Hubbard model in a harmonic trap 

Ultracold fermionic atoms on optical lattices represent nowadays the most 
direct experimental realization of the Hubbard model, since atoms interact 
only via a contact potential. In addition, in the experiments all parameters can 
be controlled with an unprecedented precision. Since the particles are trapped 
by a confined potential, it is possible to drive the MMIT by changing the total 
filling of the trap, the on-site repulsive interaction, or varying the curvature 
of the confining potential. The Hamiltonian in this case can be written as 

H = ~ l E ( c tc i+lff + H.c.) +UY, mimi + V2 E x i n ^ ( 2 ) 

i,(T i i(r 

where V2 is the curvature of the harmonic confining potential, and Xi measures 
the position of the site i (xi = ia with a the lattice constant). The number 
of lattice sites is N, and is selected so that all the fermions are confined in 
the trap. We denote the total number of fermions as Nf and consider equal 
number of fermions with spins up and down (N^ = Nf± = Nf/2). In our 
simulations, we used the zero-temperature projector QMC method described 
in detail in Refs. [21,22,23]. 

Results for the evolution of the local density (n, = (n^ + n^)), as a function 
of the total number of the confined particles, are shown in Fig. 2. For the 
lowest fillings, so that n < 1 at every site, the density shows a profile with the 
shape of an inverted parabola, similar to that obtained in the non-interacting 
case [24], and hence, such a situation should correspond to a metallic phase. 
Increasing the number of fermions a plateau with n = 1 appears in the middle 
of the trap, surrounded by a region with n < 1 (metallic). Since in the periodic 
case, a Mott insulator appears at n — 1, it is natural to identify the plateau 
with such a phase. The Mott-insulating domain in the center of the trap 
increases its size when more particles are added, but at a certain filling this 
becomes energetically unfavorable and a new metallic phase with n > 1 starts 
to develop in the center of the system. Upon adding more fermions, this new 
metallic phase widens spatially and the Mott-insulating domains surrounding 
it are pushed to the borders. Depending on the on-site repulsion strength, 
they can disappear and a complete metallic phase can appear in the system. 
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Fig. 2. Evolution of the local density in a parabolic confining potential as a function 
of the position in the trap and increasing total number of fermions. The parameters 
involved are N = 100, U = 6t and V^a 2 = 0.006i. The positions are measured in 
units of the lattice constant a. 

Finally, a "band insulator" (i.e., n = 2) forms in the middle of the trap for the 
highest fillings. Due to the full occupancy of the sites, it will widen spatially 
and pushes the other phases present in the system to the edges of the trap 
when more fermions are added. 

Alternatively to the case shown in Fig. 2, where the transitions between the 
phases were driven by an increment of the filling in the trap, Mott-insulating 
regions can be obtained by increasing the ratio U/t. Fig. 3(a) shows the evo- 
lution of the density profiles in a trapped system when this ratio is increased 
from U/t — 2 to U/t — 8. It can be seen that for small values of U/t (U/t = 2) 
there is only a metallic phase present in the trap. As the value of U/t (U/t = 4) 
is increased, a Mott-insulating phase tries to develop at n — 1 while a metallic 
phase with n > 1 is present in the center of the system. As the on-site repulsion 
is increased even further (U/t = 6, 8), a Mott-insulating domain appears in 
the middle of the trap suppressing the metallic phase that was present there. 
In Fig. 3(b) we show the variance of the density for the profiles in Fig. 3(a) 
(from top to bottom, the values presented are for U/t = 2, 4, 6, 8). As ex- 
pected, the variance decreases in both the metallic and Mott-insulating phases 
when the on-site repulsion is increased. When the Mott-insulating plateau is 
formed in the density profile, a plateau with constant variance appears in the 
variance profile with a value that will vanish only in the limit U/t — ► oo. The 
dependence of the variance in the Mott insulator vs U/t is depicted in Fig. 
3(d) up to U/t = 20 (five times the band- width). In this figure it is possible 
to see that after a fast decrease up to around U/t = 8, the variance reduces 
slowly [proportionally to (U/t)~ 2 ] when increasing U. As shown in Fig. 3(b), 
whenever a Mott-insulating domain is formed in the trap, the value of the 
variance on it is exactly the same as the one for the Mott-insulating phase 
in the periodic system for the same value of U/t (horizontal dashed lines). 
This would support the validity of the commonly used local density (Thomas- 
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Fermi) approximation [25]. However, the insets in Fig. 3(b), show that this is 
not necessarily the case, since for U/t = 4, the value of the variance in the 
Mott-insulating phase of the periodic system is still not reached in the trap, 
although the density reaches the value n — 1. Therefore, in contrast to the 
periodic Mott-insulating region is not only determined by the filling. 

In the cases of U/t = 6 [inset in Fig. 3(b) for a closer look] and U/t = 8, the 
value of the variance in the periodic system is reached and then one can say 
that Mott-insulating phases are formed there. 



n 




Fig. 3. Profiles for a trap with V20 2 = 0.00254 and Nj = 70, the on-site repulsions 
are U/t =2 (A), 4 (y), 6 (O), and 8 (O)- (a) Local density, (b) variance of the 
local density, (c) local compressibility k as defined in Eq. (3), (d) variance of the 
density in the Mott-insulating state (when U/t > 0) vs U/t. The dashed lines in 
(b) are the values of the variance in the n = 1 periodic system for U/t = 2, 4, 6, 8 
(from top to bottom). 



Although the variance gives a first indication for the formation of a local Mott 
insulator, an ambiguity is still present, since there are metallic regions with 
densities very close to n = and n = 2, where the variance can have even 
smaller values than in the Mott-insulating phases. Therefore, an unambiguous 
quantity is still needed to characterize the Mott-insulating regions. We pro- 
posed a new local compressibility as a local-order parameter to characterize 
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the Mott-insulating regions[14,15], that is defined as 




Xi,i+. 

\j\<W) 



3 ' 



(3) 



where 



Xij = {riifij) - (m) (rij) 



(4) 



is the density-density correlation function and £{U) ~ b£(U), with £(U) the 
correlation length of Xi,j m the unconfined system at half-filling for the given 
value of U. As a consequence of the charge gap opened in the Mott-insulating 
phase at half filling in the periodic system, the density-density correlations 

x 

decay exponentially [X{x) ^ ex P e(C7) ] enabling £(U) to be determined. The 
factor b is chosen within a range where k 1 becomes qualitatively insensitive to 
its precise value (b ~ 10) [15]. Physically, the local compressibility defined here 
gives a measure of the change in the local density due to a constant shift of 
the potential over a finite range but over distances larger than the correlation 
length in the unconfined system. 

In Fig. 3(c), we show the profiles of the local compressibility for the same 
parameters as Figs. 3(a) and (b) (we did not include the profile of the local 
compressibility for U = 2t because for that value of U we obtain that I is 
bigger than the system size). In Fig. 3(c), it can be seen that the local com- 
pressibility only vanishes in the Mott-insulating domains. For U = 4t, it can 
be seen that in the region with n ~ 1 the local compressibility, although small, 
does not vanish. This is compatible with the fact that the variance is not equal 
to the value in the periodic system there, so that although there is a shoul- 
der in the density profile, this region is not a Mott insulator. Therefore, the 
local compressibility defined here serves as a genuine local order parameter to 
characterize the insulating regions that always coexist with metallic phases. 



4 Local quantum criticality and universality 

In the previous section we have characterized the local Mott-insulating and 
metallic phases quantitatively. We study in this section the regions where 
the system goes from one phase to another. Criticality can arise, despite the 
microscopic spatial size of the transition region 2 , due to the extension in 
imaginary time that reaches the thermodynamic limit at T = 0, very much 
like the case of the single impurity Kondo problem [26], where long-range 
interactions in imaginary time appear for the local degree of freedom as a result 
of the interaction with the rest of the system. An intriguing future question, for 

2 Recent experiments leading to a MMIT [8,9] considered a systems with linear 
dimensions ~ 100a, i.e. still in a microscopic range. 
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both theory and experiment, will be the role of spatial dimension in the critical 
behavior of systems in the thermodynamic limit. In the periodic case, the FC- 
MMIT in ID and 2D belong to different universality classes, with dynamical 
exponents z = 2 (ID) and z = 4 (2D) [20]. This route to the transition in the 
periodic case, is the one relevant for the trapped systems where the density 
changes on entering in the local Mott insulator, and could lead to different local 
quantum critical behavior between traps with different dimensionalities. In the 
bosonic case the FC-MMIT exhibits the same dynamical exponent (z = 2) 
in all dimensions [27,28], which could also lead to a different local quantum 
critical behavior for trapped bosons as compared with trapped fermions. 
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Fig. 4. The local compressibility k vs 5 = 1 — n at 5 — > for (A) Nf = 70, U = 8t 
and V 2 a 2 = 0.0025*; (y) N f = 70, U = 6t and V 2 a 2 = 0.0025*; (Q) N f = 72, 
U = 6t and a quartic potential with V^a 4 = 1.0 x 10 -6 £; (O) unconfined periodic 
system with U = 6t. The straight line displays a power-law behavior w = 0.72. 
Inset: Dependence of the critical exponent w on At 2 . 

Figure 4 shows the local compressibility vs 5 = 1 — n for 5 — * in a double 
logarithmic plot. A power law k 1 ~ 5 m is obtained, with to < 1, such that 
a divergence results in its derivative with respect to n, showing that critical 
fluctuations are present in this region. Since the QMC simulation is affected 
by systematic errors due to discretization in imaginary time, it is important to 
consider the limit Ar — >• in determining the critical exponent. The inset in 
Fig. 4 shows such an extrapolation leading to w ~ 0.68 — 0.78. At this point we 
should remark that the presence of the harmonic potential allows the deter- 
mination of the density dependence of various quantities with unprecedented 
detail on feasible system sizes as opposed to unconfined periodic systems, 
where systems with 10 3 — 10 4 sites would be necessary to allow for similar 
variations in density. In addition to the power law behavior, Fig. 4 shows that 
for 5 — > 0, the local compressibility of systems with a harmonic potential but 
different strengths of the interaction or even with a quartic confining potential, 
collapse on the same curve. Hence, universal behavior as expected for critical 
phenomena is observed also in this case. This fact is particularly important 
with regard to experiments, since it implies that the observation of criticality 
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should be possible for realistic confining potentials, and not only restricted to 
perfect harmonic ones, as usually used in theoretical calculations. However, 
Fig. 4 shows also that the unconfmed case departs from all the others. Up to 
the largest systems we simulated (600 sites), we observe an increasing slope 
rather than the power law of the confined systems. Actually, we observe that 
the exponent of the power law obtained between contiguous points in Fig. 4 
for the periodic case extrapolates to w — 1, as it is shown in Fig. 5. 

1 

0.95 
03 0.9 
0.85 
0.8 

0.003 0.006 0.009 0.012 

(1-n) 

Fig. 5. Exponent w of the local compressibility for the periodic system, obtained 
between contiguous points in Fig. 4, vs (1 — n). The extrapolation leads to w = 1 
for n—*l. 

Having shown that the local compressibility displays universality on approach- 
ing a Mott-insulating region, we consider the variance A as a function of the 
density n for various values of U and different confining potentials. Figure 6 
shows A vs n for a variety of systems, where not only the number of particles 
and the size of the system are changed, but also different forms of the confin- 
ing potential were used. Here we considered a harmonic potential, a quartic 
one, and a superposition of a harmonic, a cubic and a quartic one, such that 
even reflection symmetry across the center of the system is broken. It appears 
at first glance that the data can only be distinguished by the strength of the 
interaction U, showing that the variance is rather insensitive to the form of 
the potential. The different insets, however, show that a close examination 
leads to the conclusion that only near n = 1 and only in the situations where 
at n = 1 a Mott insulator exists, universality sets in. The inset for n around 
0.6 and U — 8t, shows that the unconfmed system has different variance from 
the others albeit very close on a raw scale. This difference is well beyond the 
error bars. Also the inset around n — 1 and for U = 4t, shows that systems 
that do not form a Mott-insulating phase in spite of reaching a density n — 1, 
have a different variance from those having a Mott insulator. Only the case 
where all systems have a Mott-insulating phase at n — 1 (U = 8t), shows 
universal behavior independent of the potential, a universality that encom- 
passes also the unconfmed systems. For the unconfmed system, the behavior 
of the variance can be examined with Bethe-Ansatz [29] in the limit 5 — > 0. 




9 
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Fig. 6. Variance A vs. n for (O) harmonic potential V 2 a 2 = 0.00252; with N = 100; 
(A) quartic potential V4C1 4 = 5 x 10~ 7 t with N = 150; (y) harmonic potential 
V 2 a 2 = 0.016* + cubic ^a 3 = 1.6 x 10" 4 t + quartic F 4 a 4 = 1.92 x W~ 5 t with 
N = 50; and (full line) unconfined periodic potential with N = 102 sites. The 
curves correspond from top to bottom to U/t = 0,2,4,8. For a discussion of the 
insets, see text. 

In this limit and to leading order in 5, the ground state energy is given by 
[30] E (6)/N — E (5 = 0)/N oc 5, such that the double occupancy, which can 
be obtained as the derivative of the ground-state energy with respect to U, 
will also converge as 6 towards its value at half-filling. Such behavior is also 
obtained in our case as shown in Fig. 7 for the same parameters of Fig. 4. 
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Fig. 7. Variance of the density vs<5 = l — nat<5^0 for (A) Nj = 70, U = 8t and 
V 2 a 2 = 0.0025t; (y) N f = 70, U = 6t and V 2 a 2 = 0.0025t; (Q) N f = 72, U = 6t 
and a quartic potential with V^a 4 = 1.0 x 10 -6 i; (O) unconfined periodic system 
with U = 6i. The straight line displays linear behavior, i.e. its slope is equal to one. 
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5 Momentum distribution function 



In most experiments with quantum gases carried out so far, the MDF, which is 
determined in time-of-flight measurements, played a central role. A prominent 
example is given by the study of the superfluid-Mott-insulator transition [8] 
in the bosonic case. As shown below, we find that the MDF is not appropriate 
to characterize the phases of the system in the fermionic case, and does not 
show any clear signature of the MMIT. 
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Fig. 8. Normalized momentum distribution function for U/t = 2 (a) 4 (b), 6 (c), 8 
(d), and N = 100, N f = 70, V 2 a? = 0.0025*. 

In Fig. 8 we show the normalized momentum distribution function for the 
density profiles shown in Fig. 3. For the trapped systems, we always nor- 
malize the MDF to be unity at k — 0. We first notice that n& for the pure 
metallic phase in the harmonic trap [Fig. 8(a)] does not display any sharp fea- 
ture corresponding to a Fermi surface, in clear contrast to the periodic case. 
The lack of a sharp feature for the Fermi surface is independent of the pres- 
ence of the interaction and is also independent of the size of the system. In the 
non-interacting case, this can be easily understood: the spatial density and the 
momentum distribution will have the same functional form because the Hamil- 
tonian is quadratic in both coordinate and momentum. When the interaction 
is present, it could be expected that the formation of local Mott-insulating 
domains generates a qualitatively and quantitatively different behavior of the 
momentum distribution, like in the periodic case where in the Mott-insulating 
phase the Fermi surface disappears. In Fig. 8(c), it can be seen that there is 
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no qualitative change of the MDF when the Mott-insulating phase is present 
in the middle of the trap. Quantitatively in this case is similar to the pure 
metallic cases Fig. 8(a),(b). Quantitative changes in ny. appear only when the 
on-site repulsion goes to the strong-coupling regime, but this is long after the 
Mott-insulating phase has appeared in the system. 

At this point one might think that in order to study the MMIT using the 
MDF, it is necessary to avoid the inhomogeneous trapping potential and use 
instead a kind of magnetic box with infinitely high potential on the bound- 
aries. However, in that case one of the most important achievements of the 
inhomogeneous system is lost, i.e., the possibility of creating Mott-insulating 
phases for a continuous range of fillings. In the perfect magnetic box, the 
Mott-insulating phase would only be possible at half filling, which would be 
extremely difficult (if possible at all) to adjust experimentally. The other possi- 
bility is to create traps which are almost homogenous in the middle and which 
have an appreciable trapping potential only close to the boundaries. This can 
be studied theoretically by considering traps with higher powers of the trap- 
ping potentials. As shown below, already non-interacting systems make clear 
that a sharp Fermi edge is missing in confined systems. 




-500 -250 250 500 jt/4 7t/2 371/4 1 
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Fig. 9. Exact results for Nf = 840 noninteracting trapped fermions in a lattice with 
1000 sites and a confining potential Vioa 10 = 7 x W~ 27 t. Density profile (a) and the 
normalized momentum distribution function (b): the continuous line corresponds 
to (a), the dashed line is the result when an alternating potential V a = O.lt is 
superposed on the system, and the dotted line corresponds to V a = 0.5i. 

In Fig. 9(a) we show the density profile of a system with 1000 sites, Nf = 840, 
and a trapping potential of the form Vioxj with V^roa 10 = 7 x 10 -27 t. It can 
be seen that the density is almost flat all over the trap with a density of the 
order of one particle per site. Only a small part of the system at the borders 
has the variation of the density required for the particles to be trapped. In 
Fig. 9(b) (continuous line), we show the corresponding normalized momentum 
distribution. It can be seen that a kind of Fermi surface develops in the system 
but for smaller values of k, is always smooth and its value starts decreasing 
at k = 0. In order to see how changes when an incompressible region 
appears in the system, we introduced an additional alternating potential, so 
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that in this case the new Hamiltonian has the form 

H = -t£( c U + i CT + H.c.) + V 10 5>!° n la + V a J2(-lfn ia , (5) 

i,a itr ia 

where V a is the strength of the alternating potential. For the parameters pre- 
sented in Fig. 9(a), we obtain that a small value of V a (V a = O.lt) generates a 
band insulator in the trap, which extends over the region with n ~ 1 (when 
V a = 0). However, the formation of this band insulator is barely reflected in 
rik, as can be seen in Fig. 9(b) (dashed line). Only when the value of V a is 
increased and the system departs from the phase transition [V a = 0.5t, dotted 
line in Fig. 9(c)], does a quantitatively appreciable change in appear. 



6 Phase Diagram 

Finally we consider the phase diagram of the system. As shown in Figs. 2 and 
3, unlike the exponents, phase boundaries seem to be rather sensitive to the 
choice of potential, number of particles and strength of the interaction. As in 
the unconfmed case, we would expect to be able to relate systems with dif- 
ferent number of particles and/or sizes by their density. Given the harmonic 
potential, a characteristic length (in units of the lattice constant) is given by 
C = (V2 /t) , such that a characteristic density can be defined. Figure 10 
shows that the characteristic density p = iV//C is a meaningful quantity to 
characterize the phase diagram. There, the phase diagrams for two systems 
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Fig. 10. Phase diagram for a system with N = 100, V^a 2 = 0.006t (y) and N = 150, 
V2C1 2 = 0.002t (O) sites. The phases are explained in the text. 

with different sizes (N = 100 and iV = 150) and different strength of the har- 
monic potential (V^a 2 = 0.006t and V^a 2 = 0.002t respectively) are depicted 
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showing that such a scaling allows to compare systems with different sizes, dif- 
ferent number of particles, and different strength of the potential. With this, 
it is possible to relate the results of numerical simulations to much larger ex- 
perimental systems. The different phases obtained are: A pure metal without 
insulating regions (A), a Mott-insulator at the center of the trap (B), a metal- 
lic intrusion at the center of a Mott-insulator (C), a "band insulator" (i.e. 
with n — 2) at the center of the trap surrounded by a metal (D), and finally a 
"band insulator" surrounded by a metal, surrounded by a Mott-insulator with 
the outermost region being again a metal (E). Two features are remarkable 
here. The first one is that on varying the filling of the trap, a reentrant behav- 
ior is observed for the phase A. The density profile shows a shoulder as can 
be seen in Fig. 2 before reaching the plateau with n = 2 but, as shown by the 
inset of Fig. 3 for U = 4 around n = 1, it is possible to go through a region 
with n = 1 without reaching the value of the variance that corresponds to a 
Mott-insulator. The second intriguing feature is that the boundary between 
the regions A and B remains at the same value of the characteristic density 
for all values of U that could be simulated. The latter feature is consistent 
with the fact that for U — > oo density properties of two-component fermionic 
systems become identical to the ones of single species fermions with the same 
total filling. In this case, the characteristic density for the formation of the 
insulating state with n = 1 in the middle of the trap was obtained in Ref. [31] 
p ~ 2.6 — 2.7. This is the same value observed in Fig. 10 for the formation of 
the Mott insulator. A comparison between density profiles for systems with 
a Mott insulator in the middle of the trap for U/t = 8 and U/t = oo is pre- 
sented in Fig. 11(a). No big differences are observed between both cases. This 
is clearly in contrast to the comparison of the variance profiles [Fig. 11(b)] 
where in the Mott plateau the variance has changed abruptly. 
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7 Conclusions 



We have studied ground state properties of two-species fermions confined on 
ID optical lattices, reviewing and enhancing the analysis in Refs. [14,15]. We 
have defined a local-order parameter, that we called local compressibility, 
which characterizes the local insulating regions in an unambiguous way. It 
vanishes in the insulating phases and is finite in the metallic ones. This local 
compressibility gives a measure of the local change in density due to a constant 
shift of the potential over finite distances larger than the correlation length of 
the density-density correlation function for the Mott-insulating phase in the 
unconfined system. We found that the local compressibility exhibits critical 
behavior (a power-law decay k ~ |1 — n| ro ) on entering the Mott-insulating 
regions. Due to the microscopic nature of the phases, spatial correlations ap- 
pear not to contribute to the critical behavior discussed there. This is a new 
form of metal-Mott-insulator transition not observed so far in simple periodic 
systems, and that might be realized in fermionic gases trapped on optical lat- 
tices. The exponent of the power-law (zu) was obtained to be independent of 
the confining potential and/or strength of the interaction, excluding, however, 
the unconfined case. Universal behavior was also observed for the variance 
A ~ |1 — to | when to — > 1. In this case, the observed behavior is shared by the 
unconfined model. 

As opposed to periodic systems, where the appearance of the gap in the insu- 
lating phase can be seen by the disappearance of the Fermi surface in the mo- 
mentum distribution function, we found that due to the presence of a harmonic 
confining potential such a feature is much less evident. In a non- interacting 
case, we have shown that although increasing the power of the confining po- 
tential sharpens the features connected with the Fermi edge, it remains quali- 
tatively different from the homogeneous case. Hence, in fermionic systems the 
momentum distribution function does not seem to be the appropriate quantity 
to detect experimentally the formation of Mott-insulating domains. 

Finally, we determined a generic form for the phase diagram that allows to 
compare systems with different values of all the parameters involved in the 
model. It can be used to predict the phases that will be present in future 
experimental results. The phase diagram also reveals interesting features such 
as reentrant behavior in some phases when parameters are changed, and phase 
boundaries with linear forms. Results obtained for finite values of the on-site 
repulsion were contrasted with the ones for infinite U. We observe that once 
a Mott plateau has appeared in the middle of the trap, density profiles do 
not change much by increasing the on-site repulsion. This is in contrast to 
the variance of the density, which reduces continuously to zero in the Mott- 
insulating phase. 
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